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ABSTRACT 


The problem of roll, sway and yaw motions of surface ships is considered. A 
mathematical model is developed which consists of the nonlinear maneuvering equations 
and incorporates cross coupling between sway force, yaw moment and the roll angle induced 
during a steady turn. The hydrodynamic derivatives and coefficients of a typical container 
ship were used as the base-line study model. The coupled system of nonlinear algebraic 
equations is formulated and solved to predict the steady state roll angle, sway velocity and 
turning rate as a function of the rudder angle. The results are then compared to that of the 
decoupled systems currently employed. A local perturbation is implemented in the vicinity 
of the above steady states to investigate dynamic stability of motion. Sensitivity analysis with 
respect to important design parameters such as speed loss during turning, approach speed, 
transverse metacentric height and trim is performed. Results demonstrate the significance 
of the coupling between roll, sway and yaw and the need to incorporate similar studies in 


the ship design and analysis process. 
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I. INTRODUCTION 


A. GENERAL 

Any study of the motion of a ship through water requires the use of models. 
The model could be a mathematical representation which attempts to describe the 
situation symbolically, searching for a quantitative solution, or a physical model 
which strives to replicate or simulate the actual conditions for the purpose of 
collecting empirical data. The two can be combined such that they complement each 
other, or most often, physical models are used to test the validity and accuracy of 
mathematical models. However, both models fail to completely describe the physical 
phenomenon. The physical model falls short by the inability to simultaneously satisfy 
all the parameters describing the problem. The mathematical model suffers by 
departures from the true state through assumptions made during the development 
of the equations [Ref. 1]. Notwithstanding these shortcomings, the use of 
models in predicting a ships response to its environment, i.e., wind and waves, and 
its own propulsive and control systems is essential in the design process to ensure 
safety of the ship, performance of its mission, its survivability in extreme conditions 
and its efficiency during normal operating conditions. [Ref. 1] 

In general, the examination of ship motions is separated into: 1) steering and 


maneuverability (calm water), and 2) seakeeping (waves, current), with motion 


stability (absence of external excitation) and control (external excitation) some of the 
primary concerns. (Ref. 1] 
The basis of the mathematical model is rigid body dynamics, hence Newton’s 


laws of motion: [Ref. 2] 


~ ———_—__—__——_—__» 
F = £ x (linear momentum) 

(1.1) 
M = < x (angular momentum) 


The ship is free to rotate (roll, pitch, yaw) and translate (surge, sway, heave) in all 
six degrees of freedom and the forces and moments acting on the vehicle are 
comprised of: (Ref. 2] 

1. Hydrodynamic forces and moments on the bare hull, appendages, rudder and 

propeller. 

2. Inertial reaction forces and moments. 

3. Wind, waves and currents. 

4. External forces. 
Moreover, these forces and moments are dependent upon properties of the rigid 
body (e.g., geometry, mass, center of gravity), its orientation (i.e., to inertial frame 
of reference and to the body fixed frame of reference), its dynamics (velocities, 
accelerations, propeller speed, rudder deflection and deflection rate), and the fluid 


properties (e.g., density, viscosity, pressure, energy). (Ref. 1] 


At this point, the mathematical model has rapidly developed into a very 
complex and cumbersome system and requires simplification in order to be of 
practical use. Choice of the body fixed coordinate system to take advantage of any 
symmetries is one simplification used almost exclusively. Eda, et al, [Ref. 3] 
found in one study of bulbous bow type ships that the slight asymmetry of the 
underwater hull form increased significantly in roll such that this simplification was 
not valid. Another reduction in complexity often used is to separate the lateral 
motions (sway, yaw, roll) from the longitudinal motions (surge, pitch, heave) 
(Refs. 4, 5, 6]. Later studies [Refs. 3, 7, 8, 9, 10] account for the loss 
of speed during a turn by including the surge equation with the lateral motions for 
maneuvering studies. Seakeeping studies [Refs. 8, 11, 12, 13] heretofore have 
separated out the roll equations using a one degree of freedom system only, coupling 
effects were neglected as very small or incorporated as forcing terms. Rutgersson and 
Ottosson [Ref. 8] superimposed the maneuvering model motions with the seakeeping 
motions. 

As can be ascertained from this discussion, the forces, moments, velocities and 
accelerations acting upon the vessel are highly nonlinear, involving complex coupling 
of terms. Linear theory has been successful in predicting and analyzing directionally 
Stable ships with controls fixed for small perturbations [Ref. 14]. In addition, the 
linear hydrodynamic forces and moments have been reduced to semi-empirical 


equations based on ship design characteristics such as block coefficient, length, 


breadth and draft [Refs. 2, 15]. However, there has been no completely analytical 
procedure which predicts the nonlinear hydrodynamic forces and moments [Ref. 2]. 

It is necessary therefore, to combine the mathematical model with the physical 
model such that the physical model is used to measure these nonlinear forces and 
moments and incorporate the results in the mathematical model. "The functions that 
May approximate the hydrodynamic forces and moments can be expressed formally 
by Taylor expansion around the state of equilibrium with respect to the quantities 
affecting the different motions, such as the axial speed uw, turning rate r, rudder angle 
6, etc. The polynomial coefficients thus describing the hydrodynamic forces and 
motions can be determined from captive model tests at different kinds of basic 
motions [Ref. 8]. " A complete treatment of the captive model test procedure can be 
found in Principals of Naval Architecture, Vol. III [Ref. 2], and will not be discussed 
here. The Taylor expansion and subsequent model development is elaborated in the 
following Chapters. 

Captive model tests for the lateral motions are more difficult than for the 
longitudinal motions due to the large sizes required to adequately represent the 
physical model. Consequently, there exists little data for the third order coefficients 
necessary to produce semi-empirical equations for the nonlinear hydrodynamic forces 
and moments arising from these lateral motions similar to those which exist for the 
first order terms. One study by Inoue, et al, [Ref. 15] has produced formulations for 
limited nonlinear terms in yaw and sway based on drift angle and turning rate. The 


yaw moment coefficients so produced are not in good agreement with the model test 


data under all conditions of loading and motions induced. Another study done by 
Son and Nomoto [Ref. 9], utilizes model test data for the third order hydrodynamic 
coefficients of a SR-108 container ship. In the absence of such semi-empirical 
formulations for the forces and moments in the equations of motion and access to 
model test facilities, the aforementioned data, together with vehicle design 
characteristics, can be used in a mathematical model to ascertain particular 


information about the coupling of yaw, sway and roll motions. 


B. SCOPE OF THIS THESIS 

This research is intended as a bridge between the static roll restoring moment, 
controls fixed analysis of the coupled lateral motion problem and a fully coupled 
treatment of a ships roll response to maneuvering in a seaway. As discussed 
previously, current studies in the area of maneuvering and steering range from three 
degree of freedom models, to four degree of freedom models which account for the 
speed loss during a turn. These studies, however, use the equilibrium point in the 
Taylor expansion of the representative hydrodynamic forces and moments, as straight 
upright, zero rudder deflection, such that v.>=1r,>=¢)=6,=0 and all perturbations and 
nonlinearity effects are evaluated based on this single nominal point. This research 
has as its basic premise that the nominal point for the Taylor expansion changes 
significantly with rudder angle. As a result of this change, the stability characteristics 
of the ship may also change appreciably, a fact which needs to be established for 


successful prediction of ship maneuvering response in a seaway. 


Therefore, the objective is to determine the significance of the coupling of roll 
into sway and yaw motions as a function of rudder angle. In analyzing this coupled 
problem, the effects of roll into sway/yaw and conversely, sway/yaw into roll were 
investigated. 

Development of the mathematical model for the coupled yaw, sway, and roll 
equations of motion along with models for the rudder forces and moments, the 
righting moment, and the hydrodynamic forces and moments is presented in Chapter 
II. The final form of the model being reduced to a set of four nonlinear, coupled, 
first order differential equations. The additional equation is derived from kinematics 
to relate the roll angle and the roll rate. 

In order to utilize Taylor expansion for formulation of the functions which 
approximate the hydrodynamic forces and moments, the state of equilibrium for the 
motion parameters, v, r, and @, is established. This is accomplished through steady 
state analysis. In Chapter III, this analysis is first done for the linear case to obtain 
an initial approximation for the case of near zero rudder deflection (6~0). This 
linear steady state solution is incorporated into the nonlinear steady state analysis 
in an iterative process which produced the equilibrium values for yaw rate, sway 
velocity, and roll angle as a function of the rudder deflection in a steady turn. As the 
rudder angle is incrementally increased, the steady state values from the previous 
rudder angle are used as the initial approximations to initiate the iterations for the 


current rudder deflection. 


Once the steady state solutions have been computed, their stability properties 
are established by local perturbation of the coupled equations of motion in the 
neighborhood of each nominal point. This procedure, described in Chapter IV, yields 
a generalized eigenvalue problem. The solutions to the eigenvalue problem provide 
information as to the stability properties of the nominal points and thus characterizes 
the stability of the coupled sway, yaw and roll motions. Demonstration of the 
coupling effect on the lateral motions, is shown utilizing comparison to two 
independent systems: 

1. A sway/yaw model represented by two nonlinear equations decoupled from 
roll in that the hydrodynamic forces and moments are functions of sway 
velocity (v) and yaw rate (r) only. 

2. A roll model represented by the nonlinear roll equation and the kinematic 
roll rate equation decoupled from sway and yaw in that the hydrodynamic 
forces and moments are functions of roll angle (@) and roll rate (p) only. 

The directional and roll stability characteristics in a steady turn are established 
by the degree to which the real part of the complex eigenvalues is negative and the 
magnitude of the roll damping ratio. 

Chapter V relates the sensitivity and stability of design parameters such as 
metacentric height and trim to the coupling between the lateral motions. Results are 
presented in terms of rudder angle with a qualitative summary at the end of this 
chapter. 


Finally, Chapter VI summarizes the conclusions and presents recommendations 


for further research and utilization of the mathematical model herein developed. 


Il. MATHEMATICAL MODEL 


A. EQUATIONS OF MOTION 


For ship motions, the full nonlinear equations of motion in a body fixed axis 


system are shown in Equations (2.1). These equations reflect an assumption of 


symmetry along the longitudinal axis (1.e., y,=0). The orthogonal, right-hand axes for 


the body fixed system is shown in Figure 2.1. The parameters are defined in Table 


II. 


>< 


Surge: m[u -x(q? +r?) +Z(q+pr) -rv+wq] 


Tl 
< 


Sway: m[v+x.(pq+r)+z,(qr-p) tur-wp] 


Heave: m[w+x,(rp-4)-z,(p?+q*) +pv-qu] = Z 


[p+(,-1,)qr-1,(@q +r) -mz,(v+ur-wp) 
) = K-AGZ() 


Roll 


,q-(.,-1,)rp-1,,(r?-p?) +mlzg(i+qw-1v)- 
X_(w+pv-qu)]=M 


Pitch: 


Yaw: [7 +Q,-L)Pq +I _(qr-p)+mx(v+ur-wp) = N 


(2.1a) 


(2.1b) 


(2.1c) 


(2.1d) 


(2.1e) 


(2.1f) 





Figure 2.1: Coordinate System [Ref. 11] 
B. SIMPLIFYING ASSUMPTIONS 


For the purpose of this investigation, a mathematical model having three 
degrees of freedom vice the six degrees of freedom depicted in Equations (2.1) was 
used. The following simplifying assumptions were made: 

1. The rotational velocity and acceleration about the y-axis are zero. (q = 0 
and q=0) 


2. The translational velocity and acceleration in the z direction are zero. (Ww 
= 0 and w=0) 


3. The vertical heave and pitch motions are decoupled from the horizontal plane 
motions. 


4. The product of inertia I, is very small and can be neglected. 


5. The surge equation is substituted by an algebraic equation which is a function 
of u, V, and 6. 


6. The longitudinal center of gravity, (LCG) and the longitudinal center of 
buoyancy, (LCB) are at midship. 


7. The vertical center of gravity, (VCG) is on the centerline. 

8. The only important forces and moments acting on the ship induced by the 
rudder are those due to rudder deflection. Forces and moments due to 6 and 
5 are negligible. [Ref. 2] 


TABLE I. EQUATIONS OF MOTION PARAMETERS 










| Parameter | Deseription 
Distance along the principal axes 


u,V,W Translational velocity components of ship relative to fluid | 

along body axes | 
Rotational velocity components of ship relative to inertial 
P.q.r reference system along body axes 


- . = a 
* 


=e 


a | 
i 


Se 


Hydrodynamic force components along body axes 
Hydrodynamic moment components along body axes 


w yaw angle: bow to starboard positive 
8 pitch angle: bow up positive 
¢@ roll angle: starboard down positive 


Mass of ship 


Xo,¥q.ZcG Coordinates of the center of gravity in the body axis 
system 





| A | Displacement weight of ship 
GZ(@) 





Rudder angle in radians 







Initial velocity of ship 


Vv 
| cess 2 heel Mass density of sea water 


Ship length between perpendiculars (LBP) 
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Applying these assumptions to Equations (2.1), the three degree of freedom 


equations of motion for the model used are given in Equations (2.2). 


Sway: m[v+x f-zptur] = Y (2.2a) 
Roll: [p-mz,(vtur) = K- AGZ(4) (2.2b) 
Yaw: I7+mx,(vtur) = N (ZZ) 


Equations (2.2) are nondimensionalized for ease of working between model test 
data and actual ship test data using the relationships shown in Table II. 
To demonstrate the use of these nondimensionalizing terms, substitution of the 


appropriate values from Table II into Equation (2.2a) yields: 


] V? |. ivan mor V ] 
Eee als ie ne 
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TABLE II. NONDIMENSIONAL PARAMETER RELATIONSHIPS 


os 








AGZ($)’ = 


AGZ(o)| — 


Joey 
7? 


12 


simplifying by factoring out (V’/L): 


(3 ot2v') m!t9!sxg¢/-2gp"su'r’) = (3pt2v") ¥ 


and Equation (2.2a) becomes Equation (2.3) in nondimensional form. 


m'[ v! + xor/ - zop! + u'r’) = Y’ (2.3) 


Since Equation (2.3) and Equation (2.2a) are of the same form, the prime notation 
will be dropped and all equations will be considered as represented in 


nondimensional form unless otherwise indicated. 


C. FORCE AND MOMENT REPRESENTATION 
Using the method of Abkowitz and Strom-Tejsen [Ref. 2], the sway force, roll 
moment and yaw moment can be expressed as Equations (2.4). 
Sway Force: Y =f, (u,V,7,U,V,F,,0,0) (2.4a) 
Roll Moment. K = f,(u,v,r,u,¥,7,0,6,0,8) (2.4b) 


Yaw Moment: N = f,(u,v,r,,v,7,0,6,) (2.4¢c) 


A third order Taylor expansion of f,, f,, and f,; can be expressed as Equations 
(2.5). 
fr = Yyt+7p+Y we +Y vr2+¥ r+¥ 7? +¥ vert 


Y,O+Y,O+Y,.V°O+Y vb +Y ro? +Y, arb +¥(3)+Y, 


(2.5a) 
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f, = K,v+K,v3+K,,V°rtK y+K, 734K, vr?+K,o+K,o+K,o+ (2.5b) 


K mg +K wah +K eg +K(5)+Ky-A GZ(o) 
fy = NytN +N, vetN, vr7+N tN s+N_p2+N vert (2.50) 
.SC 
NiO+N O+N VO +N Vb7+N, r+ rh? +N(5) +o 
The coefficients in Equations (2.5) are the hydrodynamic coefficients and 
represent the partial derivative with respect to the subscripted variable. For example, 
Y, means hs and Y,,, means acl #Y . The terms Yj, Ky, and N, are the sway force, 
ov OV\ or? 
roll moment and yaw moment induced by the propeller respectively. In a similar 
manner Y(5), K(d) and N(6) represent the force and moments induced by the 
rudder. The righting moment due to the static stability of the ship is represented by 


the AGZ(¢) term. 


D. RUDDER FORCE AND MOMENT REPRESENTATION 
The expressions used to determine the rudder force and moments were taken 
from Son and Nomoto [Ref. 9] and are presented in Equations (2.6) and (2.7). The 


parameters are defined in Table III. 


Y(5) = - (1 + a,)F,cos(d) (2.6a) 
K(5) = (1 + a,)z,F,,cos(5) (2.6b) 
N(8) = - (Xp + @y%,)Fycos(8) (2.6c) 
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TABLE III. RUDDER FORCE AND MOMENT PARAMETERS 


Parameter Description 
Rudder to hull interaction coefficient 


“4 Normal force action on the rudder 
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z coordinate of point on which rudder force Y;, acts 
x coordinate of point on which rudder force Y; acts 
x coordinate of point on which normal force F,, acts 
Rudder aspect ratio 

Rudder area 


R 


R Effective rudder inflow velocity 


=) 


Hl Effective rudder inflow angle 


VR Components of rudder effective inflow velocity 
constant in Equation (2.7c) 

Effective propeller inflow velocity 

constant in Equation (2.7c) 


Thrust coefficient 


= 


Advance coefficient 
Number of propeller revolutions per second 


Propeller diameter 


< 


Effective propeller wake fraction 


cS 


constant in Equation (2.7e) 


x coordinate of propeller position 


Cenc. propeller flow rectification coefficients 


7 


CoC. Cy rudder wake coefficients 


flow rectification coefficients 
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6.13A \ {4 
F, = —j| V,* sin(a 
us ery Fa x Sin(st,) 
Vp = Up se Vp’ 
8kK, 
Up = Une 1 + ; 
nJ 
7; - “2a 
nD 


u, = u[(l-w,)+t(v+x 7)? +c vtec, ri] 


Ve TA estoy ey av 
v 
a, = & + an 
Up 


RIGHTING MOMENT REPRESENTATION 


(2.7a) 


(2.7) 


(2.7c) 


(2.7) 


(2.7e) 


(2.7£) 


(2.78) 


For the initial part of the righting arm curve and for wall sided ships, Equation 


(2.8) gives a good approximation for the function GZ(@). [Ref. 16] 


GZ() = GMsin(d) + =~ tan*(6) sin() 


GM = the transverse metacentric height 
where BM = the transverse metacentric radius 
@ = the roll angle 


(2.8) 


Equation (2.8) expresses the fact that for most ships of fairly rectangular 
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midship section, the GZ(¢) curve exhibits the typical characteristics of a hardening 


spring restoring moment. The spring constant is an increasing function of the roll 
angle @. To see this, we use the Taylor series expansion up to third order for the 


trigonometric functions: 


sin(p) 


tl 
& 
| 


tan(>) 


ll 
& 
lL 


Therefore, if terms higher than third order are neglected: 


Z 
tan?()sin() = (¢ 30 [-26" : (6-20) a 


and substituting into Equation (2.8) yields: 
_ Gules 243) . R743 - Gy lay — 1lapy)g3 


Now for most ships, 3 BM/>GM, Which means that the leading cubic coefficient 


ie the GZ(¢) curve is positive and the curve concaves up as would a hardening 


spring. 


F. STATE-SPACE REPRESENTATION 
Combining Equations (2.2) and (2.4) and rearranging the terms, the equations 


of motion can be expressed in vector form: 
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a F(z, 5) 
(2.9) 


where z is the state vector [vr p $]’. 
Equation (2.9) is more conveniently expressed in matrix notation as is shown in 


Equations (2.10). The last matrix equation is derived from the kinematic roll rate. 


(mx,-N,)  (I,-N,) O -N, 


v Y 
(m-¥,) (mxg-¥,) -mzg -¥%y| [Fi _ Us (2.10a) 
-MZ¢ 0 (,-K,) -K,| |P R 
| oO 0 0 1 d 
= Nv+N_v?+N_vr2+(N_- +N _7>+N_v2r+N,+ 
ty v wv wT ( r MX gU)r rr wr a 9? (2.10b) 
Nig? OtN yg VO° +N, gh O+N, gr? +N(5) +No 
= Yv+Y y34+¥ vr2+(¥ -mu)r+Y_r°+Y. v?r+Y,6+ 
ts VARY KV, Mr +Y, Ft ¥ Vv THY (2.10c) 
YgY O+¥ g4Vb°+¥ gr O+Y yyrd + ¥(5)+¥, 
= Ky+K_v°+K._ y2r+(K.+m K_r°?+K_vr?+K,o+ 
Tr v vw wr" ( r Zgu)r+ rf Ry VE s? (2.10d) 


Kg¥ O+K yg VO° +K, gO +K 44rd’ +K(5)+Ky-A GZ(9) 
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III. STEADY STATE REPRESENTATION 


A. LINEAR STEADY STATE SOLUTIONS 

With the above third order Taylor expansion, the equations of motion 
represented by Equations (2.10) are highly nonlinear, the solution of which depends 
upon the initial and equilibrium conditions. The steady state condition requires the 


acceleration terms y ¢ 5, @ to be equal to zero. Thus Equation (2.10a) reduces to 


Equation (3.1) for steady state: 


“| = [0] (3.1) 


where fy, f;, and fp are defined in Equations (2.10). 
The numerical solution to Equation (3.1) is an iterative process which depends 
upon a fairly accurate initial estimate of the root for convergence. This starting point 


can be determined by assuming an equilibrium state of v)>=r1,.>=p,=¢)=0 and a very 


small initial rudder deflection (6,<1°). Application of these assumptions to Equation 


(3.1) yields the linear, steady-state, yaw-sway-roll equations of motion in Equation 


(3.2). 
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N, (N,-mxgu) —N, - -N(3)-N, 


Y, (¥,-mu) im r| = |-¥(8)-Y, (3.2) 
K, (K,+mz gu) (K,-4GM)| |? -K(5)-K, 


Equation (3.2) is in the familiar matrix form Az = B, where z is the state vector 
with p=0. 

Various computer subroutines are commercially available which solve this 
simple system of equations for the state vector values, given the other parameters. 
The values for the state vector obtained in this manner are then used as the initial 


estimate in solving the nonlinear system of equations. 


B. NONLINEAR STEADY STATE SOLUTION 

Once the initial estimate for equilibrium is obtained, the system of nonlinear 
equations in Equation (3.1) can be solved numerically for the state vector nominal 
points for each rudder angle, 6. Here again, commercial computer subroutines are 
available. ‘The Levenberg-Marquardt algorithm with analytic Jacobian was chosen as 
the solution method [Ref. 17]. Since p=0 in Equation (3.1), the system of 
nonlinear equations is reduced to three and the Jacobian was determined by 


Equations (3.3). 
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(3.3a) 


S 
Y/Y 2B eB 
21% VB g|s 
ZL ZlLs ZS 


a = N,+3N,,V° +N, 77+2N, wr t2N, vb Nyy? SO (3.3b) 


oe =2N,,,vr+(N,-mxgu)+3N,_7°+N,,v°+2N,_ rh +N 44h? a (3.3c) 


2 = Ni +N,4¥7 +2N avd +N, gr? +2N, rd + — (3.3d) 

2 = Y+3Y v?+V 7?+2¥ wr+2Y meV tY sag? Sate 

2 =2Y,,vr+(Y,-mu)+3Y_77+Y,v?+2Y ro+y,., ?+ ee) (3.3f) 

Ys = Yet Y cv" +2Y 44¥0 spp Melle) ER rd + —— en) (3.3g) 
op 

os = K,+3K,,v?+K,7?+2K,,vr+2K,,.V0+K,.,0°+ a“ (3.3h) 


—*-2K,,vr+(K,+mzgl)+3K, J? +K,,¥?+2K rg ro +K yo? SO) el) 
r 


SR AK, +K, v+2K .v+K 72 +2K ar -KAGLD) —(3,3)) 


op 
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Nb) = Fr) a 
where (3) = Fv) = a - 0  (3.3k) 
K(8) = Fyvy) aK). 
op 
and 


— =A GMcos(6)BMtan*(6)( ——— 3e05(0)] oar 


Since F,, F,, and F; are complicated functions of v and r (Equations (2.6) and 


(2.7)), the following approximation method was used for evaluation: 


Lg I ie ca (3.4) 
ae eee x, - 0.99x, 

Fach term in the Jacobian is determined by evaluation at the previous nominal 
point with the initial nominal point determined from linear solutions as described in 
Section A above. 

The result of the numerical solution to the nonlinear equations (Equations 


(3.1)) is the steady state values for v, r, and @ for each rudder angle under steady 


turning conditions and for a given GM and propeller speed. 


C. TYPICAL RESULTS 
Utilizing the solution method depicted in Sections A and B above, PROGRAM 


COUPLED listed in Appendix A, was developed to predict the steady state roll 
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angle, sway velocity and turning rate of a high speed container ship as a function of 
rudder angle. 

The design characteristics of a SR-108 container ship are presented in Table 
IV. The hydrodynamic derivatives and coefficients used in the computer simulation 
are listed in Table V. 

For a base-line model, the hydrodynamic derivatives and coefficients were held 
constant in addition to the Froude number (hence the propeller revolutions), the 
transverse metacentric height, and speed loss ratio. Therefore, for the figures 


presented in this section: 


Froude number Fn = 0.3 
Metacentric height GM = 0.3 m 
Speed loss ratio a = 0.6 


Figures 3.1 through 3.3 show the variation of sway velocity, yaw rate, and roll 
angle as a function of rudder angle during a steady turn. It can be seen that all 
Steady state variables v, r and @ are highly nonlinear functions in 6 which is 
attributed to both the effect of the nonlinear terms in our equations of motion and 
the fact that the ship speed is reduced during the turn. Had a constant speed, linear 


model been utilized, v, r and @ would appear as straight lines versus 6. 
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TABLE IV. PRINCIPAL DIMENSIONS OF SR-108 CONTAINER SHIP [Ref. 9] 


Hull Length B. P. 
Breadth 
Draft Fore 
Aft 
Mean 
Displacement Volume 
Height from keel to 
transverse metacenter 
Height from keel to 
center of buoyancy 
Block coefficient 
Prismatic coefficient 
Waterplane area 
coefficient 
Midship section 
coefficient 
LCB from forward 
Perpendicular 


Bilge keel 
Length 
Depth 


Rudder Area 
Height 
Aspect Ratio 
Area Ratio 


| Propeller 
Diameter 
Pitch Ratio 
Expanded Area Ratio 
Boss Ratio 
Number of blades 
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33.0376 
7.7583 


3.000000 
0.435000 
0.137100 
0.154300 
0.145700 
0.106860 


0.17810 


0.07912 
0.55900 
0.58000 
0.68600 
0.518 L 
0.240 L 


0.009709 

0.133000 

1.821900 
1/45.0 





TABLE V. HYDRODYNAMIC DERIVATIVES AND COEFFICIENTS [Ref. 9] 


A.) Hull Only 


0.0214000 -0.0424000 -0.00055800 
ace Prat ae 
Tae Pas ee 


NY 0.0093250 Nie -0.0038592 Ke -0.00024300 : 
thd 


Td 
Y ibd -0.0013560 0.0024195 0.00003569 


B.) Propeller and Rudder 


(1, 
|cx | 070 | & | oo | | 05% 
| % | 0033 | cy | 0156 | Ke | 0.s2704sst 





0.000 C= -0.275 V 0.088 v>0 
| 0.193 v<0 | 


79.10 Fn=0.2 
118.64 Fn=0.3 
158.19 Fn=0.4 





2S 


v. SWAY VELOCITY 





SET -20 10 0 10 20 30 


RUDDER ANGLE (deg 
Figure 3.1: Sway Velocity Under Steady Turning 


r: YAW RATE 





RUDDER ANGLE (deg 
Figure 3.2: Turning Rate Under Steady Turning 
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Figure 3.3: Roll Angle Under Steady Turning 
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IV. STABILITY ANALYSIS 


A. COUPLED STEERING AND ROLL EQUATIONS OF MOTION 
Once the nominal points are determined for each rudder angle, Equation 


(2.10a) can be reduced to state-space representation (Equations (3.4)) for stability 


analysis. 
BZ = Az 
v v 
: ' (4.1a) 
where g = |. : z= 
Pp Pp 
>| b 
(mx,-N,) (,-N,) 0 “Ny 
_ (m-Y,) (mx,-Y,) -mZ> -Y; (4.1b) 
“MZ> 0 (1,-K,) “Ky 
0 0 0 1 | 
hy hy Hy Hy 
Ov or dp od 
Gs Gs As Ys 
Fee ie: (4.1c) 
Gr Gr Gr An 
Ov or dp op 
op op op op 
Ov or dp oad 


Matrix B can be thought of as the generalized mass matrix for the problem 
while matrix A is the Jacobian matrix of f(z,6) evaluated at each nominal point. For 
stability analysis purposes, the variables v, r, p, and ¢@ are understood to represent 
small deviations of the actual variables from their respective steady state values as 
computed in the previous chapter. The rationale behind Equations (4.1) is the 
introduction of small local perturbations in v, r, p, and @ superimposed on their 
steady turning values. Lyapunov’s linearization theorem establishes that stability 
properties of a nominal point for a nonlinear system can, in general, be deduced 
from the stability properties of the corresponding linearized system. 


The terms in the matrix A are as defined in Equations (3.3) with the following 


additions: 
Op Op Op 
e _ P _ ® _ 4g 
Ov or op 
gy i 
Op 


Stability can be determined from the solution to the eigenvalue problem: 


ABz = AZ 
(4.2) 
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If all eigenvalues, A, have negative real parts, the nominal point is asymptotically 
stable. If at least one eigenvalue is positive, the nominal point is unstable. The 
degree of stability of the nominal point is therefore related to the real part of the 
eigenvalues. The more negative they are, the faster the exponential convergence of 
solutions in time to the nominal point. The imaginary part of the eigenvalues 
characterizes the frequency content of the oscillatory behavior of the system 
response. 

A subroutine was included in the computer model to solve this eigenvalue 


problem as a function of rudder angle. 


B. DECOUPLED STEERING EQUATIONS 
In order to distinguish the characteristic contributions of the sway and yaw 
motions from those of roll, the steering equations (Equations (4.3)) were decoupled 


from the maneuvering equations (Equations (4.1)), resulting in Equations (4.3). 











ies el BD (4.3a) 
(m-Y,) (mx,-Y,)| [|r - 
his = YyytY, vi +¥,_ vr? +(¥ -mu)r+Y_ 3+ v?r+¥(8) (4.3b) 
Fys = NytN ve +N, vr?+(N -mxcuyr+N_ oP +N, v?r+N(8) (4.3c) 


This decoupling is based on the assumption that all cross coupling coefficients 


between roll and sway/yaw are zero, which is the usual approximation made for 
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surface ships where roll is studied independently from sway/yaw. One of the 
motivations for this study is to evaluate the degree of accuracy of such an 
approximation. 

The additional subscript, S, distinguishes the decoupled steering forces and 
moments from those of the fully coupled system. Following the same format as 


before, the state-space formulation for stability analysis is given in Equations (4.4). 











Byis A 52s 
my : (4.4a) 
where Z. = ‘ ee os ; 
I-N 
B. (mxg-N,) (,-N,) (4.4b) 
(m-Y.) (mx,-Y,) 
as 
roe aL (4.4c) 
Gss Ass 
OV or 
ys = N+3N v7+N_7?+2N y+ ONG) (4.4d) 
av v VVV vur ov 
ts _on _vr+(N Mg) * BN r24N_y2+ N06) (4.4) 
Or or 
me = ¥, + 3Y, vy? +¥ 7242¥ vrs 9) (4.4f) 
“2 = DY, V+, MU) +BY nF? #¥ pg ¥?+ OP (4.4g) 
: 
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The approximation method of Equation (3.4) was used to determine rudder 
force and moment contributions in the Jacobian matrix A. 

The stability attributable to yaw and sway as a function of rudder angle is thus 
determined by the solution of the decoupled steering eigenvalue problem in 


Equations (4.5). 


sy (4.5) 


C. DECOUPLED ROLL EQUATIONS 
By decoupling the steering equations, the roll equations become uncoupled 


producing the following relationship: 




















UL -K, : K, 3 = (4.6a) 
0 1 | {¢ P 
fer = Kyo + K(5) - AGZ() (4.6b) 


The additional subscript, R, refers to the decoupled roll forces and moments 
and GZ(¢) is as given in Equation (2.8). Equations (4.7) are the state-space stability 


analysis formulations. 
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(4.7a) 
Beet 
where eae epics, 
(I,-K,) -K, 
Bp = , (4.7b) 
Sen Are (4.7¢) 
4. |? © 
“| @p 
ap a 
Tee gy, Tre _ x HAGA) (4.74) 
Op * Ab Op 
z = 1 : “ = Q 
p 
(4.7e) 


The expression — is given in Equation (3.3). 


Therefore, the stability attributable to roll as a function of rudder angle is 


attained by solution of the decoupled roll eigenvalue problem of Equations (4.8). 


(4.8) 
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D. TYPICAL RESULTS 

Continuing with the base-line model described in Chapter III Section C, the 
roots of Equations (4.2), (4.5) and (4.8) were obtained using PROGRAM 
COUPLED (Appendix A). Figure 4.1 is a plot of the real roots of the fully coupled 
maneuvering equations as a function of rudder angle. As can be seen, there are three 
distinct components which were classified as due to roll or due to steering (sway/yaw) 
by decoupling the steering and roll equations as discussed in Sections B and C above. 
Figures 4.2 through 4.4 show the three components of the fully coupled system 
matched with their respective decoupled solutions of Equations (4.5) for steering and 
Equations (4.8) for roll. 

In Figure 4.2, the decoupled rol] real roots are constant for each rudder angle, 
whereas the coupled roots vary with an increasing stability to about five degrees of 
rudder. The upper steering component shown in Figure 4.3 indicates very little 
deviation in the coupled and decoupled roots, however, the lower steering 
component in Figure 4.4 shows a deviation in the two roots as the rudder angle is 
increased. The steering eigenvalues become increasingly more negative as the rudder 
angle is increased. This indicates that the steering system is dynamically more stable 
for a non-zero rudder angle, and as a result, a steering system that is designed to be 
stable for straight line motions will be even more stable for motions along curved 
reference paths. 

The imaginary roots for the decoupled equations are zero, indicating that all 


imaginary roots obtained in the solution of the coupled equations are due to roll 
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only. Figure 4.5 shows that for rudder angles less than 15°, the coupled model 
produces larger values than the decoupled model and for rudder angles greater than 
15°, the reverse is true. Comparing Figures 4.2 and 4.5, we can see that the imaginary 
part of roll response, or the roll frequency of oscillation, can be predicted fairly 
accurately from the decoupled roll equations. This is not the case for the real part 
where decoupling the equations results in a severe underestimation. 

Figure 4.6 is a root locus plot of the roll eigenvalues for both coupled and 
uncoupled models. As can be seen, the fully coupled model is more stable than the 
decoupled model. 

The roll damping caused by roll is shown in Figure 4.7. The decoupled roll 
analysis shows less damping for each rudder angle than does the fully coupled model. 
Furthermore, the damping increases significantly for large rudder angles. It can be 
seen that the actual roll damping is significantly larger than what is predicted from 
the decoupled model, which may underestimate the actual damping ratio by a factor 


of two to three. 
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Figure 4.2: Coupled Model and Decoupled Model Roll Eigenvalues 
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Figure 4.3: Coupled Model and Decoupled Model Upper Steering Eigenvalues 
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Figure 4.4; Coupled Model and Decoupled Model 
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Figure 4.7: Coupled Model and Decoupled Model Damping Ratio 
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V. SENSITIVITY OF MATHEMATICAL MODEL 
TO VARIATION IN DESIGN PARAMETERS 


A. BACKGROUND 

Since actual ships rarely operate under the precise conditions assumed in the 
development of this mathematical model, changes in certain variables to measure the 
response sensitivity is prudent, as for any design study. The following sensitivity 
analysis has this aim in mind as well as to identify any inexplicable response as an 
area for further study. It must be remembered that certain assumptions were made 
at the outset which might need refinement as the design spiral continues. For 
instance, the roll damping term, Kj, is assumed to be constant throughout the 
simulation, although experimental evidence [Ref. 18] suggests that roll damping 
coefficients are forward speed and motion amplitude dependent. By isolating a 
parameter, problem areas are enhanced so that a particular parameter can be the 
focus for analysis of the validity of the assumptions made or the modelling method 
used. 

If, on the other hand, one or more of the dependent se do not show 
appreciable sensitivity, further simplifying assumptions could possibly be made when 
analyzing the more difficult problems of maneuvering in regular and irregular seas. 

In the following sections, sensitivity to speed loss during a turn, Froude number 


of approach spped, metacentric height, operating turn, and righting moment is 
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evaluated. For brevity, a limited number of plots are presented in each case. 
Appendix B contains additional plots for metacentric height, Froude number, and 


speed loss. 


B. SPEED LOSS DURING A TURN AS THE DESIGN PARAMETER 
The sensitivity of the model to changes in speed loss during the turn was 
investigated by holding transverse metacentric height (GM) and Froude number 


constant: 


In determining the speed loss under steady turning conditions, it was assumed 
that the water depth was sufficient, i.e. > 110 meters, to neglect this effect [Refs. 2, 
3]. However, it has been determined by several sources: Davidson (1944), Shiba 
(1960), Strom-Tejsen (1965), Eda and Crane (1965), among others [Ref. 2], that 
speed loss in steady turning is a function of the hull configuration (including rudder) 
and the turning diameter. Various other sources, [Refs. 3, 9, 10] consider rudder 
angle as an additional parameter in determining speed loss. In an effort to combine 
these effects, Figures 25 and 178 and Tables 29(a) and 29(b) from Principles of Naval 
Architecture [Ref. 2], Figure 11 of Eda and Falls [Ref. 3], and Figure 10 of Son and 
Nomoto [Ref. 9] were used to develop the following simple relationship to represent 


speed loss: 


41 


= l-a*6, 


© 


wa 
V, 
V = steady turning speed 


V, = approach speed 
where 
a = speed loss ratio 


5, = rudder angle in radians 


Variation of the speed loss parameter ranges from no speed loss (a#=0) to maximum 
speed loss (a=0.73) during a turn with fixed controls. Results of these effects on 
sway velocity, yaw rate and roll angle may be viewed in Figures 5.1 through 5.3. The 
influence of speed loss on r is negligible, indicating that in all cases considered, the 
ship steady turning rate performance is not affected. On the contrary, v and ¢ show 
a relatively significant modification due to speed loss. The sign reversal of the roll 
angle plot (Figure 5.3) for large rudder angles is rather uncommon, but it is, 
nevertheless, possible as a result of the nonlinearities that are present in the 
equations of motion for this model. 

The characteristic roots as a function of rudder angle were determined to 
establish stability trends. The real roots indicate three distinct components as 
previously discussed in Chapter IV, Section C. Both the upper and lower steering 
components demonstrate a significant increase in stability with rudder angle but little 
variance in stability due to speed loss. The roll component, however, indicates 
greater stability as speed loss decreases. We discovered earlier, that the non-zero 
values in the root locus were due solely to the roll contribution owing to the 


imaginary parts of the steering elements. The damping ratio of the roll contribution 
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imaginary parts of the steering elements. The damping ratio of the roll contribution 
is related to both real and imaginary parts of the eigenvalues. Figures 5.4 and 5.5, 
for the fully coupled and uncoupled models respectively, show this as a function of 
rudder angle. Both models exhibit qualitatively the same trend, with the uncoupled 


model showing a significant underestimation of the true damping ratio value. 


C. FROUDE NUMBER AS THE DESIGN PARAMETER 
In testing the model sensitivity to Froude number changes, the transverse 


metacentric height and speed loss ratio were kept constant: 


GM = 03m 
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Figure 5.4: Effect of Speed Loss on Coupled Model Damping Ratio 
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Froude numbers used to calculate the forward velocity prior to initiating the turn 
aswell as propeller revolutions were assumed to remain constant during the turn. For 
each Froude number and rudder angle, the sway velocity, yaw rate and roll angle 
were determined using the model program. Figures 5.6 through 5.8 at the end of this 
section depict graphic results. Note that v and r are not readily affected by change 
in Froude number, whereas ¢ Is greatly changed, increasing dramatically with Froude 
number. 

Using these values, the stability of the coupled response was determined. The 
real roots of the coupled maneuvering equations for the three Froude numbers are 
indistinguishable when plotted as in Figure 4.1, as are the roots for the lower 
steering components (Figure 4.2). The upper steering component exhibits nonlinear 
behavior through five degrees of rudder for all Froude numbers with increasing 
stability as both Froude number and rudder angle increase. The roll component was 
nonlinear throughout the range of rudder angles tried. Furthermore, as the Froude 
number increased, so did the range of rudder angles for which stability increased. 

The coupled model roll damping ratio remains relatively insensitive to Froude 
number (Figure 5.9), while the uncoupled model suggests a wide variation (Figure 
5.10). This demonstrates the incorrect conclusions which can be reached when using 
linear/uncoupled models that neglect the actual coupling of the lateral sway/yaw 


motions back into roll. 
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Figure 5.10: Effect of Fn on Decoupled Model Damping Ratio 


D. TRANSVERSE METACENTRIC HEIGHT AS THE DESIGN PARAMETER 
For transverse metacentric height as the design parameter in the model, Froude 


number and speed loss ratio were held constant: 


Fn = 0.23 
0.6 


6 


Here again, the sway velocity, yaw rate and roll angle for each GM and rudder 
angle were computed using the model program. The results are represented in 
Figures 5.11 through 5.13. As indicated for the case of Froude number being the 


design parameter, v and r are affected very little by even an order of magnitude 
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change in GM. On the other hand, the increase in GM not only causes a large 
decrease in the magnitude of the roll angle, but in the range as well. 

In determining the influence of GM on the stability characteristics, it was found 
that the real roots plotted as a function of rudder angle showed no distinction 
between GM’s. When the components were separated out, only the upper steering 
and roll constituents exhibited significant variation. The upper steering response 
reflected a somewhat incongruous reaction, that is, as the GM increases, the stability 
due to steering decreases slightly. The roll component demonstrated a marked 
increase in stability with increase in GM. This result is evident in the root locus 
diagram of Figure 5.14. 

Conversely, the damping ratio did not change significantly, as shown in Figure 
5.15. Computation of the damping ratio from the decoupled model showed a much 


wider deviation among different GM’s. 


E. TRIM AS THE DESIGN PARAMETER 
For the investigation of a vessel in the trimmed condition, the transverse 


metacentric height, Froude number and speed loss ratio were held constant: 


GM = 03m 
Fn = 0.3 
a = 0.6 
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The base-line model chosen for this study, and therefore the hydrodynamic 
derivatives and coefficients derived from model tests, were for a design trim of 0.5 
meters by the stern and a mean draft of 8.5 meters for the actual ship (refer to Table 
IV for model parameters). It has been shown [Ref. 11] through linear theory that 
these derivatives and coefficients change for the trimmed condition by the following 


semi-empirical relations: [Ref. 15] 




















: F fa 
Y(t) = YO)1 a 
Na) = N,0 5 
Y(t) = YO + i 
N,(t) = a “| 
— 
. e 
where _ + for trim by stern 
-— - for trim by bow 


d_ = design mean draft 


Y(O) = Y, at design condition 


The application of these relationships are somewhat questionable, but in the 
absence of model tests, they will suffice to indicate trends for further study. 
A further modification in the model was made to the zz parameter (Table III) 


or the z coordinate of the point on which the rudder force acts: 
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d 
z(t) = z,(0) + tale 


Figures 5.17 through 5.19 depict the effects of ra = +0.2 which equates to a 


significant change in operating condition. The sone of v and r to the previous 
parameters was practically negligible with the exception of the speed loss ratio. This 
is not the case for trim where turning rate is improved for trim by the bow and 
decreased for trim by the stern. It should be noted that this is contrary to linear 
theory results. For sway velocity at small angles of rudder the variation is small, but 
for large rudder angles the trim by stern curve deviates considerably from the other 
two conditions. The difference is more pronounced in the roll angle curves (Figure 
5.19) where trim by the bow generates higher angles during the turn. On the other 
hand, trim by the stern results in smaller roll angles with a roll angle sign reversal 
beyond a certain rudder angle. This is due to the fact that for large rudder angles, 
the hydrodynamic moment exerted by the rudder force dominates over the fluid 
forces on the hull. 

The root locus diagram (Figure 5.20) for the roll stability analysis is revealing 
in that very different trends can be observed. Trim by the stern is considerably more 
stable than by the bow for the entire range of heel angle, although at larger rudder 
angles, the stability is less than at the design condition. This is possibly qualitatively 
explained by noticing in Figures 5.21 and 5.22, that the damping ratio for trim by the 


stern increases sharply precisely at the inflection point of the roll angle as shown in 


5 


Figure 5.19. Here again, it must be pointed out that the damping ratio’s variance 
between the coupled model and the decoupled model 1s at least a factor of two. It 
can be seen from both coupled and decoupled models that at the inflection point of 
the roll angle, the damping ratio reaches its peak value. This maximum value of the 
damping ratio exists for the other trim conditions as well, but its actual location is 


shifted outside the range of rudder angle variations considered here. 


F. RIGHTING ARM CURVE AS THE DESIGN PARAMETER 
The sensitivity of the model to changes in the righting arm as a function of the 
roll angle was tested by holding the transverse metacentric height, Froude number 


and speed loss constant: 


GM = 10m 
Fn = 0.3 
a = 0.6 


In general, the righting arm curve, GZ(@), is taken to be a very simplified 
version of Equation (2.8), namely GZ(¢)=GM*¢@. This linear righting arm gives 
good approximations for small angles of heel and is derived much the same as was 
done for Equation (2.8) in Chapter II, Section E. For the purpose of this 
investigation, the simplified GZ(¢) was left in its nonlinear form, 
GZ(¢)=GM*sin(@). Figure 5.23 illustrates the two GZ(@) curves used. Note that 
GZ] (representing Equation (2.8)) develops a much larger righting arm at increased 


angles of heel than does GZ2 (representing GM*sin(@)). 
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Figure 5.21: Effect of Trim on Coupled Model Damping Ratio 
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The effects of these different righting arm curves on the model for v, r, and @, 
are shown in Figures 5.24 through 5.26. Sway velocity and yaw rate appear to be 
unaffected by this parameter change in contrast to roll angle, which shows relatively 
high sensitivity, as expected. 

Figures 5.27 (coupled model) and 5.28 (decoupled model) indicate the same 
trend in damping ratio sensitivity as the previous parameters, specifically, the 


decoupled model underestimates ¢ by a factor of two. 


G. SUMMARY OF RESULTS 
A qualitative summary of results is presented in Table VI below. All 


comparisons are to the base-line configuration of:’ 


GM = 03m 
Fn = 0Q3 
a = 0.6 


' For GZ sensitivity GM = 1.0 m. 
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Notation in the table refers to the largest noted deviations from the model 


according to the following subjective scale: 


Symbol Standard 


Seats — 
Negligible < 1% 

O 1% - 5% 

Moderate 5% - 20% 
High > 20% 


Mixed => large change 
over the range 





Additionally, A, B and C refer to rudder angle range for the noted deviations such 


that: 


0 <A < 10 
10 < B < 20 
20 < C < 30 
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TABLE VI. QUALITATIVE SENSITIVITY SUMMARY 


Parameter 
V 
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This reference table shows clearly that the sway velocity and turning rate are 
relatively insensitive to Froude number, metacentric height and righting moment, 
while being quite sensitive to the speed loss and trim parameters. Roll and damping 


ratio show the greatest response for all the parameters considered. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

This thesis developed a mathematical model for the coupled sway, yaw and roll 
equations of motion with rudder deflection as the independent variable. Additionally, 
various components of the hydrodynamic forces and moments, e.g., rudder to hull 
interaction, were modelled using existing research formulation. Deviation of this 
study from current solution methods was in the manner of evaluating the Taylor 
expansion for the steady state Jacobian matrix. The equilibrium position for each 
rudder angle was found to deviate significantly from the commonly used zero state, 
even for rudder angles of less than 10°. In fact, the roll angle and sway velocity 
nominal points are nearly maximum at 10° of rudder. 

Comparison of stability results for the coupled lateral motion model to those 
of models representing the steering equations of motion decoupled from the roll 
equations of motion showed: 

1. The directional stability characteristics signified by the real part of the 
complex eigenvalues are not affected significantly by the coupling of roll into 
sway and yaw. This indicates that sway and yaw can be treated separately 
from roll when studying directional stability. 

2. Roll stability represented by the magnitude of the damping ratio is 
considerably affected by the coupling of sway and yaw into roll. The fact that 


the damping ratio is underestimated by a factor of two by the uncoupled 
treatment of sway and yaw, indicates ship design by this means is a significant 
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over-design that bears more attention, especially with the economic 
constraints of today’s Navy. 


3. The natural frequency is not significantly affected by coupling. This indicates 
that while maneuvering in a seaway Studies require the coupling of the lateral 
motions, the induced roll from the static righting moment can be studied from 
decoupled systems. 

4. Steering stability increases with increasing rudder deflection. This may 
indicate that if a ship is directionally unstable for straight line motion, it 
becomes directionally stable during a turn through the coupling of the lateral 
motions. 

The sensitivity of certain design parameters to the coupled system was 
investigated with the following determinations: 

1. The parameters tested were speed loss during a turn, transverse metacentric 
height, Froude number of initial forward speed, trim and righting moment. Of 
these parameters, trim proved to be the most affected. This indicates that 
during all phases of design, operating trim must be considered in addition to 


the design trim. 


2. Sway velocity and turning rate were relatively insensitive for all parameters 
except for trim and speed loss. 


3. Roll angle and damping ratio of all the parameters were significantly affected. 
4. The coupled model damping ratio showed significantly lower sensitivity than 
the decoupled mode] for smal] rudder angles in all parameters. 
B. RECOMMENDATIONS 
The results produced in this thesis were as a consequence of experimentally 
reported nonlinear hydrodynamic coefficients and design parameters of a SR-108 


container ship [Ref. 10]. Additional studies using similar data are required. 
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Utilizing the method of determining equilibrium points for each rudder 
position, design a control law for an autopilot and determine the effects of the 
coupled lateral motions on antopilot dynamics. 

The sharp change in roll angle (¢) for small values of rudder angle (6) suggests 
that this ship would be an excellent candidate for rudder roll stabilizing techniques. 

Finally, as indicated in Chapter I, Section B, use the coupled sway, yaw and roll 
equations of motion and build a mathematical model which treats wave effects in 
regular and irregular seas. Once this model is developed, analyze maneuvering in 
waves from the standpoint of biased roll oscillations. Biased roll oscillations are 


meant as roll about non-zero heel angle. 
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APPENDIX A 
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PROGRAM COUPLED 


This program solves the coupled yaw-sway-roll ship maneuvering 
equations. In addition, the yaw-sway steering eduagians are de- 
coupled from the roll equations and solved in order to compare 
results with the fully eat results. The program is broken 
into several subroutines, both original and standard commercial 
programs. 


+ &€ & & He & HH He HR H 


KAKA KEKKKKKKKKKEKKKKKEKEKKEKEEKKKEKEEEKKEKKEKEKEKKKKKEKKEKKKEKEKEKKEKKEKKKKKKK 
KEKKKKKKKEEKKKKKEKKKEKEKEKEKKKEKKEKKKEKKKKEKEEKEKKKEKEKEKEKEKEKEKEEKEKEKEKEEEEKEEEKEEKKKEKEKEEKEKKEK 


* 


* 
* 
* 


The INCLUDE statement is used to pass all the variables between 
the many subroutines. 


INCLUDE ’HYDCOMN. for’ 

COMMON /WORKSP/RWKSP 

REAL RWKSP(5000) 

REAL*8 VDEL(300) ,V(300) ,R(300) ,RPHI (300) 
EXTERNAL FCN,LSJAC,DNEQNJ 


Open the input and output files 


OPEN(10, FILE 
OPEN(15, FILE 
OPEN(20, FILE 
OPEN(25, FILE 
OPEN(30, FILE 
OPEN(26, FILE 


*COEF.INP’ ,status=’old’ ) 
*VAR.OUT’ ,status=’new’ ) 
*NLIN.DAT’ ,status=’ new’ ) 
’EIGVAL1.DAT’ ,status=’ new’ ) 
*EIGVAL2.DAT’ ,status=’ new’ ) 
*EIGVALL.DAT’ ,status=’ new’ ) 
OPEN(27, FILE =’EIGVALp.DAT’ ,status=’ new’ ) 
OPEN(35, FILE =’ZETA.DAT’ ,status=’ new’ ) 
OPEN(36, FILE =’OMEGA.DAT’ ,status=’ new’ ) 
OPEN(37, FILE =’ZROLL.DAT’ ,status=’new’ ) 
OPEN(38, FILE =’RUDDER.DAT’ ,status=’ new’ ) 
PI=4.0D0*DATAN(1.0DO) 

ERRREL=0.0001 

ITMAX=200 


This subroutine reads the hydrodynamic derivative and coefficient 
values from COEF.INP and writes pertinent values in VAR.OUT. 


CALL INPUT 
DEL=DEL* (PI/180.0D0) 


This subroutine sets up the linear EOM’s in matrix form: 
Ax=B 


CALL LINEAR 
This 1s an IMSL MATH/LIBRARY Subroutine [Ref. 19] 
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* 
* 
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Purpose: Solve a real general system of linear equations with 
iterative refinement 


Arguments: 
N -Number of equations. (Input) 
A -N by N matrix containing the coefficients of the linear 


system. (Input) 

LDA -Leading dimension of A. (Input) 

B -Vector of length N containing the right-hand side of the 
linear system. (Input) 

IPATH -Path indicator, if = 1 the system A*x=B is solved 

X -Vector of length N containing the solution to the linear 
system. (Output) 


CALL DLSARG(N,A,3,B,1,XGUESS) 

WRITE(15,*)’ THE SOLUTION TO THE LINEAR SYSTEM OF EQUATIONS IS? 
WRITE(15,*)? ? 

WRITE(15,*) ’VO=’ ,XGUESS(1) 

WRITE(15,*) ’RO=’ , XGUESS(2) 

WRITE(15,*) ’PHIO=’ , XGUESS(3) 

WRITE(15,*)? ? 

PHIO=XGUESS(3) *(180.0D0/PI) 

WRITE(15,*)’PHIO IN DEGREES =’ ,PHIO 

WRITE(15,*)? ? 


Set the rudder limits for steady turning 


DELMIN=0.01D0 
DELMAX=30 .ODO 


This Do-Loop calculates the values of sway velocity, yaw rate and 
roll angle for each rudder angle. 


DO 100 IDEL=1, IMAX 

WRITE (*,*) IDEL, IMAX 

RDEL=DELMIN+ ( DELMAX - DELMIN) *REAL ( ( IDEL-1) / (IMAX-1)) 
DEL=RDEL* (PI/180.0D0) 

UND=1 .0DO-SLR* (DABS(DEL) ) 


This is an IMSL MATH/LIBRARY Subroutine [Ref. 17] 


Purpose: Solve a system of nonlinear equations using the 
Levenberg-Marquardt algorithm with a user-supplied 
Jacobian. 
Arguments: 
FCN - User-supplied SUBROUTINE to evaluate the system of 


equations to be solved. The usage is 

CALL FCN (X, F, N), where 

X - The point at which the functions are 
evaluated. (Input) 
X should not be changed by FCN. 

F - The computed function values at the point X. 
(Output) 

N - Length of X, F. (Input) 

FCN must be declared EXTERNAL in the calling program. 

LSJAC - User-supplied SUBROUTINE to evaluate the Jacobian at a 
point X. The usage is 
CALL LSJAC (N, X, FJAC), where 
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N - Length of X. (Input) 

X - The point at which the function is evaluated. 
(Input) 
X should not be changed by LSJAC. 

FJAC - The computed N by N Jacobian at the point X. 
(Output) 


LSJAC must be declared EXTERNAL in the calling program. 


ERRREL - Stopping criterion. (Input) 
The root is accepted if the relative error between two 
successive approximations to this root is less than 
ERRREL 

N - The number of equations to be solved and the number 
of unknowns. (Input) 

ITMAX - The maximum allowable number of iterations. (Input) 


The maximum number of calls to FCN is ITMAX*(N+1). 
Suggested value = 200. 

XGUESS - A vector of length N. (Input) 
XGUESS contains the initial estimate of the root. 


X - A vector of length N. (Output) 
X contains the best estimate of the root found by 
DNEQNJ. 


FNORM - A scalar which has the following value, 
F(1)**2+...+F(N)**2 at the point X. (Output) 


CALL DNEQNJ (FCN ,LSJAC,ERRREL ,N, ITMAX,XGUESS, X , FNORM) 
PHI=X(3)*180.0D0/PI 
WRITE (20 ,50)RDEL,X(1) ,X(2) ,PHI 

50 FORMAT(4E15.5) 


Calculate the rudder forces and moments: Equations 2.6 and 2.7 


UP=UND* (WP+TAU* ( (X(1)+XP*X (2) ) **2+CPV*X(1)+CPR*X(2) )) 
JP=(U*UP) / (NP*D) 

KT=KTA+ (KTB*UP) 

UR=UP*EPSILON* ((1.0D0+8.ODO*K*KT) / (PI*(JP**2)))**0.5 

VR= GAMMA®X (1) +CDR*X (2) +CDRRR* (X (2) **3) +CORRV*X (1) * (X(2)**2) 
VD=(UR**2+VR**2) **0,5 

ALPHAR=DEL+DATAN (VR/UR) 
FLAM=(6.13D0*LAMBDA) / (LAMBDA+2.25D0) 

FNN=FLAM* (AR/ (L**2) )*( (VD) **2) *DSIN(ALPHAR) 


YD=- (1.0D0+AH) *FNN*DCOS (DEL) 
=~ (XR+AH*XH) *FNN*DCOS (DEL) 
KD=(1.0DO+AH) *ZR*FNN*DCOS (DEL) 

WRITE (38,*)ND, YD,KD 


Rudder angle in degrees 


VDEL ( IDEL) =DEL 
V(IDEL)=X(1) 
R(IDEL) =x (2) 


Rudder angle in radians 


RPHI ( IDEL) =X(3) 
XGUESS(1)=X(1) 
XGUESS (2) =X (2) 
XGUESS (3) =X (3) 
100 CONTINUE 

* 


a 


+ + + + + F & 


This subroutine solves the eigenvalue problem |A - ake | 0 
for the fully coupled yaw-sway-roll equations, using the values 
for v, r, and phi obtained above. 


CALL EIGTH1(VDEL,V,R,RPHI) 


This subroutine solves the eigenvalue problem |As - lambda*Bs|= 0 


for the steering yaw-sway equations. 


CALL EIGTH2(VDEL,V,R) 


This subroutine solves the eigenvalue problem |Ar - lambda*Br|= 0 


for the roll equations. 


CALL EIGTH3(VDEL,Rphi) 
STOP 
END 


SUBROUTINE INPUT 
INCLUDE ’*HYDCOMN. for’ 


The numbers in the variables such as N122 stand for: 
1 implies v 
2 implies r 
3 implies phi 
4 implies delta 


READ(10,*)NO,N1,N2,N3,N111,N122,N222,N112,N113,N133 ,N223, 
& N233,KO,K1,K2,K3,K111,K122,K222,K112,K113,K133,K223, K233, 
& YO, Y1,Y2,Y3,Y111,Y122,Y222.Y112.Y113,Y133, Y223 ,Y233,L ,W,XG, 
& AH; AR LAMBDA, XR, XH, ZR,BM,WP, TAU,XP, CPV, CPR, EPSILON, K,KTA,KTB, 
& D,CDR,CDRRR,CDRRV,GAMMA 


Calculate the nondimensional mass 


M=W/ (0.5D0*(L**3) ) 

WRITE(*,*)’ENTER THE GM AND THE SPEED LOSS RATIO’ 
READ (*,*)GM, SLR 

WRITE(15,*) >GM=’ ,GM 

WRITE(15,*)? 

WRITE(15,*)’ THE SPEED LOSS RATIO IS=’,SLR 


Adjusts the z coordinate of CG based on the GM 


ZGM=DABS(0.3D0-GM) /L 

ZG=0 .00078-ZGM 

WRITE(15,*)’THE SHIP DISPLACEMENT (IN CUBIC METERS) IS W=’ ,W 
WRITE (15,*)’ 

WRITE(15,*)’ THE SHIP NONDIMENSIONAL MASS IS M=’,M 
WRITE(15,*)? ? 

XG=XG/L 

WRITE(15,*)’L=’,L 

WRITE(15,*) ’XG=’ ,XG 

WRITE(15,*) ’ZG=’ ,ZG 

WRITE(15,*)? ’ 

WRITE(*,*) ENTER THE INITIAL RUDDER ANGLE IN DEGREES’ 
READ(*,*)DEL 

WRITE(15,*)’ THE INITIAL RUDDER ANGLE DEL=’ ,DEL 
WRITE(*,*) ENTER THE FROUDE NUMBER AND THE IMAX’ 
READ(*,*)FN, IMAX 
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WRITE(15,*)’ THE FROUDE NUMBER IS FN=’ ,FN 
: Calculate the initial forward speed 


U=FN* ((9.8D0*L) **0.5D0) 

WRITE(15,*)’THE FORWARD VELOCITY IS U=’,U 

WRITE(15,*)? ’ 

WRITE(*,*)’ ENTER THE PROPELLER SPEED IN RPMS FOR THE FN’ 
READ(*,*)NP 

WRITE(15,*)’ THE PROPELLER SPEED IN RPMS IS NP=’ ,NP 


* Change rpm’s to radians per sec 
NP=NP* (PI/30.0D0) 


The weight component is nondimensionalized such that WG*GM is 
dimensionless 


+ + & & 


WG= (W*9.81D0) /(0.5D0* (L**3) * (U**2 

WRITE(15,*)’THE NONDIM COMPONENT OF DISPLACEMENT IS WG=’ ,WG 
RETURN 

END 


SUBROUTINE LINEAR 


INCLUDE ’HYDCOMN. for’ 
INTEGER IROW, ICOL 
DO 100 IROW=1,N 
B( IROW) =0.0D0 
DO 90 ICOL=1,N 
A(IROW, ICOL)=0.0D0 
90 CONTINUE 
100 CONTINUE 


* 


- Equation 3.2 
* 


A(1,1)=N1 

A(1,2)=N2- (M*XG*UND) 
A(1,3)=N3 

A(2,1)=Y1 

A(2,2)=Y2- (M*UND) 
A(2,3)=Y3 

A(3,1)=K1 
A(3,2)=K2+(M*ZG*UND) 
A(3,3)=K3- (WG*GM) 


FNL=((6.13D0*LAMBDA) / (LAMBDA+2.25D0) ) * (AR/L**2) 
N4=- (XR+(AH*XH) ) *FNL 

Y4=-(1.0D0+AH) *FNL 

K4=(1.0DO+AH) *ZR*FNL 


B(1)=-(N4*DEL) -NO 
B(2)=-(Y4*DEL) - YO 
B(3)=-(K4*DEL) -KO 
RETURN 

END 


SUBROUTINE FCN(Z,Q,NDIM) 
INCLUDE ’HYDCOMN. for’ 
INTEGER NDIM 

DIMENSION Z(NDIM) ,Q(NDIM) 
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Equation 2.7 


UP=UND* (WP+TAU* ( (Z(1)+XP*Z(2) ) **2+CPV*Z(1)+CPR*Z(2) ) ) 
JP=(U*UP) / (NP*D) 

KT=KTA+ (KTB*JP) 

UR=UP*EPSILON* ((1.0D0+8.0DO0*K*KT) / (PI* (JP**2) ))**0 

VR= GAMMA®Z (1) #CDR*Z (2) #CDRRA*(Z(2)**3)+CORAV*Z(1)* (3(2)**2) 
VD=(UR**2+VR**2) **0.5 

ALPHAR=DEL+DATAN (VR/UR) 
FLAM=(6.13D0*LAMBDA) / (LAMBDA+2.25D0) 

FNN=FLAM* (AR/ (L**2))*( (VD) **2) *DSIN(ALPHAR) 


Equation 2.6 
=- (1.0D0+AH) *FNN*DCOS (DEL) 
N4=- (XR+AH*XH) *FNN*DCOS (DEL) 
K4=(1.0D0+AH) *ZR*FNN*DCOS (DEL) 
Equation 2.10b 
C1=N1*Z(1)+N111*(Z(1) **3)+N122*Z(1)*(Z(2)**2) 
C2=N222* (Z(2) **3)+N112* (Z(1) **2) *Z(2) +N3*Z(3) 
C3=N113* (Z(1) **2) *Z(3) +N133*Z(1) *(Z(3) **2)+(N2- (M*XG*UND) ) *Z(2) 
C4=N233*Z(2)*(Z(3) **2) tN4+NO+N223* (Z(2) **2) *Z(3) 
Q(1)=C1+C02+C03+C4 
Equation 2.10c 
D1=V¥1*Z(1)+Y111%*(2Z(1) **3)+Y¥122*2Z(1)*(Z(2)**2)+(Y2- (M*UND) ) *Z(2) 
D2=V222" (2(2) **3) +¥112* (2 (1) **2) *2 (2) +¥3*2Z(3) 
D3= =VII3* (Z(1) **2)*2Z(3) +¥133*2(1)* (Z(3) **2) +¥223* (Z(2) **2)*Z(3) 
D4=Y233*Z(2)*(Z(3)**2)+Y4+Y0O 
Q(2)=D1+D2+D3+D4 
Equation 2.9 


GZ=GM*DSIN(Z(3) )+0.5*BM* (DTAN(Z(3) )) **2*DSIN(Z(3) ) 
GZ=GM*DSIN(Z(3) ) 


Equation 2.10d 
E1=K1*Z(1)+K111*(Z(1) **3)+K122*2Z(1) *(Z(2) **2) +(K2+M*ZG*UND) *Z(2) 
E2=K222* (Z(2) **3) +K112*(Z(1) **2) *Z(2) +K3*Z(3) 

E3=K113* (Z(1) **2) *Z(3)+K133*Z (1) *(Z(3) **2)+K223* (Z(2) **2) *Z(3) 
E4=K233*Z(2)*(Z(3) **2)+K4+KO- (WG*GZ) 
Q(3)=E1+E2+E3+E4 


RETURN 
END 


SUBROUTINE LSJAG(NDIM, Z,QUAC) 
INCLUDE *HYDCOMN. fo 

INTEGER NDIM 

DIMENSION Z(NDIM) ,QUAC(NDIM,NDIM) 


Equation 2.7 
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UP=UND* (WP+TAU* ( (Z(1) +XP*Z(2) ) **2+CPV*Z(1)+CPR*Z(2))) 
V9=Z(1)*0.99D0 

R9=Z(2)*0.99D0 

UPV9=UND* (WP+TAU* ( (V9+XP*Z(2) ) **2+CPV*V9+CPR*Z(2) )) 
UPR9O=UND* (WP+TAU* ((Z(1)+XP*RQ) **2+CPV*Z(1)+CPR*RQ) ) 
JP=(U*UP) / (NP*D) 

JPV9=(U*UPVQ) / (NP*D) 

JPR9O=(U*UPRQ) /(NP*D) 

KT=KTA+ (KTB*JUP) 

KTV9=KTA+ (KTB*UPV9) 

KTR9=KTA+ (KTB*UPRQ) 
UR=UP*EPSILON*((1.0D0+8.0D0*K*KT) / (PI* (UP**2) ))**0.5 
URV9=UPV9*EPSILON* ((1.0D0+8.0D0*K*KTV9) / (PI* (UPV9**2)))**0.5 
URR9=UPR9*EPSILON* ((1.0D0+8.0D0*K*KTRQ) / (PI* (UPR9**2)))**0.5 
VR=GAMMA*Z(1)+CDR*Z(2) +CDRRR* (Z(2) **3) +CDRRV*Z(1)*(Z(2) **2) 
VRV9=GAMMA*V9+CDR*Z(2) +CDRRR* (Z(2) **3) +CDRRV*V9* (Z(2) **2) 
VRR9=GAMMA*Z (1) +CDR*R9+CDRRR* (R9**3) +CDRRV*Z(1) * (RO**2) 
VD=(UR**2+VR**2) **0.5 

VDV9= (URV9**2+VRV9**2) **0.5 

VDR9=(URRO**2+VRRQO**2) **0.5 

ALPHAR=DEL+DATAN (VR/UR) 

ALPHARV=DEL+DATAN (VRV9/URV9) 

ALPHARR=DEL+DATAN (VRR9/URRQ) 
FLAM=(6.13D0*LAMBDA) / (LAMBDA+2.25D0) 

FNN=FLAM* (AR/ (L**2) )*( (VD) **2) *DSIN(ALPHAR) 

FNNVO=FLAM* (AR/(L**2) ) *((VDV9) **2) *DSIN(ALPHARV) 
FNNRQ=FLAM* (AR/ (L**2) ) *((VDRQ) **2) *DSIN (ALPHARR) 


Y4=-(1.0D0+AH) *FNN*DCOS (DEL) 
Y4V9=- (1.0D0+AH) *FNNV9*DCOS (DEL) 
Y4R9=- (1.0D0+AH) *FNNRO*DCOS (DEL) 


* Equation 3.4 
* 


PY4V=(Y4-Y4V9) /(Z(1)-V9) 
PY4R=(Y4-Y4R9) /(Z(2) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS (DEL) 
N4V9=- (XR+AH*XH) *FNNV9*DCOS (DEL) 
N4R9=- (XR+AH*XH) *FNNRO*DCOS (DEL) 
PN4V=(N4-N4V9) /(Z(1)-V9) 
PY4R=(Y4-Y4R9) /(Z(2) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS (DEL) 


* Equation 3.4 
* 


N4V9=- (XR+AH*XH) *FNNVO*DCOS(DEL) 
N4R9=- (XR+AH*XH) *FNNRO*DCOS (DEL) 
PN4V=(N4-N4V9) /(Z(1) -V9) 
PN4R=(N4-N4R9) /(Z(2) -R9) 
K4=(1.0DO+AH) *ZR*FNN*DCOS (DEL) 
K4V9=(1.0D0+AH) *ZR*FNNV9*DCOS (DEL) 
K4R9=(1.0DO0+AH) *ZR*FNNRQ*DCOS (DEL) 


+ % 


Equation 3.4 


PK4V=(K4-K4V9) /(Z(1) -V9) 
PK4R=(K4-K4R9) / (Z(2) -R9) 


* Equation 3.3b 
* 
C11=N14+3*N111*(Z(1)**2)+N122* (Z(2) **2)+2.0D0*N112*Z(1) *Z(2) 
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+ * 


C12=2.0D0*N113*Z(1) *Z(3)+N133*% (Z(3) **2) 
QJAC (1,1) =C11+012+PN4V 
Equation 3.3c 


C21=2.0D0*N122*Z(1)*Z(2)+(N2-M*XG*UND) +3 .0D0*N222* (Z(2) **2) 
C22=N112*(Z(1)**2)+2.0D0*N223*Z (2) *Z(3) +N233* (Z(3) **2) 


QUAC (1,2) =C21+C22+PN4R 
Equation 3.3d 


C31=N34+N113* (Z(1)**2)+2.0D0*N133*Z(1)*Z(3)+N223* (Z(2) **2) 
C32=2.0D0*N233*Z (2) *Z(3) 


QJAC (1,3) =C31+C032 
Equation 3.3e 


D11=Y¥14+3.0D0*Y111*(Z(1)**2)+Y¥122* (Z(2) **2) +2.0D0*Y112*Z(1)*Z(2) 
D12=2 .0D0*Y113*Z(1) *Z(3)+¥133* (Z(3) **2) 


QJUAC(2,1)=D11+D12+PY4V 
Equation 3.3f 


D21=2.0D0*Y122*Z(1)*Z(2)+(Y2-M*UND) +3 .OD0*Y222* (Z(2) **2) 
D22=2 .0D0*Y223*Z (2) *Z(3)+Y¥233* (Z(3) **2)+¥112*(Z(1)**2) 


QJAC (2,2) =D214+D22+PY4R 
Equation 3.3g 


D31=Y3+Y113* (Z(1)**2)+2.0D0*Y133*Z(1)*Z(3)+Y223*(Z(2)**2) 
D32=2.0D0*Y233*Z(2)*Z(3) 


QJAC (2 ,3) =D31+D32 
Equation 3.3h 


E11=K14+3.0D0*K111*(Z(1)**2) +K122* (Z(2) **2)+2.0D0*K112*Z(1)*Z(2) 
E12=2 .0D0*K113*Z(1) *Z(3)+K133* (Z(3) **2) 


QUAC(3,1)=E11+E12+PK4V 
Equation 3.31 


E21=2.0D0*K122*Z(1)*Z(2)+(K2+M*ZG*UND) +3 .0D0*K222* (Z(2) **2) 
E22=2 .0D0*K223*Z (2) *Z(3) +K233* (Z(3) **2)+K112*(Z(1) **2) 


QUAC (3 ,2) =E21+E22+PK4R 
Equation 3.3k 
DGZ1=GM*DCOS(Z(3) )+0.5D0*BM* (DTAN (Z ( 
DGZ2=BM* (DTAN(Z(3) )*((1.0D0/DCOS(Z(3 


DGZ=DGZ1+DGZ2 
DGZ=GM*DCOS(Z(3) ) 


)**2*DC0S(Z(3) ) 


3)) 
)))**2)*DSIN(Z(3))) 


Equation 3.3) 
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E31=K3+K113*(Z(1)**2)+2.0D0*K133*Z(1)*Z(3) +K223*(Z(2) **2) 
E32=2.0D0*K233*Z(2) *Z(3) -WG*DGZ 


QUAC (3,3) =E31+E32 


RETURN 
END 


SUBROUTINE EIGTH1(VDEL,V,R,RPHI) 


INCLUDE ’HYDCOMN. for’ 
INTEGER I,J,IERR 
REAL*8 AE(4,4) ,BE(4,4) ,BETA(4) , VDEL(300) ,V(300) ,R(300) ,RPHI (300) 
REAL*8 P(300) ,EVALR(4) ,EVALI(4) ,EVAL(4) ,ALFR(4) ,ALFI(4),Z(4,4) 
EXTERNAL DMACH,RGG 
DO 500 I=1, IMAX 
P(I)=0.0D0 

SET UP THE B MATRIX AND RESET EACH I 
UP=UND* (WP+TAU* ( (V(I)+XP*R(I)) **2+CPV*V(I)+CPR*R(I) ) ) 
V9=V(I) *0.99D0 
R9=R(1)*0.99D0 
UPV9=UND* (WP+TAU* ( (V9+XP*R(I) ) **2+CPV*V9+CPR*R(I))) 
UPR9O=UND* (WP+TAU* ( (V(1) +XP*RQ) **2+CPV*V(I)+CPR*RQ) ) 
JP=(U*UP) / (NP*D) 
JPV9=(U*UPV9) / (NP*D) 
JPR9=(U*UPRQ) / (NP*D) 
KT=KTA+(KTB*JP) 
KTV9=KTA+ (KTB*JUPV9) 
KTR9=KTA+ (KTB*JUPRQ) 
UR=UP*EPSILON*((1.0D0+8.0D0*K*KT) / (PI* (UP**2) ))**0.5 
URV9=UPV9*EPSILON* ((1.0D0+8.0D0*K*KTVQ) / (PI* (UPV9**2) ))**0.5 
URR9=UPR9*EPSILON* ((1.0D0+8.OD0*K*KTRQ) / (PI* (JPRO**2) ))**0.5 
VR=GAMMA*V (I) +CDR*R (1) +CDRRR* (R(1) **3)+CDRRV*V(I)*(R(I) **2) 
VRV9=GAMMA*V9+4CDR*R (1) +CDRRR* (R( 1) **3)+CDRRV*V9* (R(1) **2) 
VRR9=GAMMA*V (I) +CDR*R9+CDRRR* (R9**3)+CDRRV*V(I) *(R9O**2) 
VD=(UR**2+VR**2) **0.5 
VDV9=(URV9**2+VRV9**2) **0.5 
VDR9=(URR9**2+VRRQ**2) **0.5 
ALPHAR=VDEL (1) +DATAN(VR/UR) 
ALPHARV=VDEL (1) +DATAN(VRV9/URV9) 
ALPHARR=VDEL (I )+DATAN(VRR9/URRQ) 
FLAM=(6.13D0*LAMBDA) / (LAMBDA+2.25D0) 
FNN=FLAM* (AR/(L**2))*((VD) **2) *DSIN(ALPHAR) 
FNNV9=FLAM* (AR/(L**2) ) *((VDV9) **2) *DSIN(ALPHARV) 
FNNRO=FLAM* (AR/(L**2) )*((VDRQ9) **2) *DSIN(ALPHARR) 


Y4=-(1.0DO0+AH) *FNN*DCOS(VDEL(1I) ) 
Y4V9=- (1.0D0+AH) *FNNV9*DCOS(VDEL 
Y4R9=- (1.0D0+AH) *FNNR9*DCOS (VDEL 
PY4V=(Y4-Y4V9) / (V(I) -V9) 
PY4R=(Y4-Y4R9) /(R(1)-R9) 

N4=- (XR+AH*XH) *FNN*DCOS(VDEL (I) ) 
N4V9=- (XR+AH*XH) *FNNV9*DCOS(VDEL (I) ) 
N4R9=- (XR+AH*XH) *FNNR9*DCOS(VDEL (1) ) 
PN4V=(N4-N4V9) /(V(I) -V9) 
PN4R=(N4-N4R9)/(R(I) -RQ) 

K4=(1.0DO0+AH) *ZR*FNN*DCOS(VDEL (I) ) 
K4V9=(1.0DO+AH) *ZR*FNNV9*DCOS(VDEL (I) ) 
K4R9=(1.0DO+AH) *ZR*FNNR9*DCOS(VDEL (I) ) 
PK4V=(K4-K4V9) /(V(I) -V9) 


(I) ) 
(I)) 


7] 


PK4R=(K4-K4R9) /(R(I) -R9) 
Calculate the fundamental circular frequency 


WN= ( (WG*GM) / (0 .000021D0))**0.5 
IF(FN.LE.O.1)TH 

KAPPA=0 . 1D0 

ELSE IF(FN.GE.0.2) THEN 

KAPPA=0 . 2D0 

ELSE 

KAPPA=FN 

ENDIF 
KPHIDOT=-0.000021D0*KAPPA*WN 


Equation 4.1b 


BE(1,1)=0.00035245D0 
BE(1,2)=0.000875D0 
BE(1,3)=0.0D0 
BE(1,4)=-0.000213D0 
BE(2,1)=0.014969D0 
BE(2,2)=0.00035245D0 
BE(2,3)=-0.000221D0 
BE(2,4)=0.0D0 
BE(3,1)=-0.000221D0 
BE(3,2)=0.0D0 
BE(3,3)=0.000021D0 
BE(3,4)=-KPHIDOT 
BE(4,1)=0.0D0 
BE(4,2)=0.0D0 
BE(4,3)=0.0D0 
BE(4,4)=1.0D0 
Equation 4.1c SET UP THE A MATRIX 


C11=N1+3*N111*(V(I) **2)+N122* (R(I) **2)+2.0D0*N112*V(I)*A(I) 
C12=2.0D0*N113*V(1I) *RPHI(I)+N133* (RPHI(I)**2) 


AE(1,1)=C11+C12+PN4V 


C21=2.0D0*N122*V(I) *R(I)+(N2-M*XG*UND) +3 .OD0*N222* (R( 
C22=N112*(V(I) **2)+2.0D0*N223*R (I) *RPHI(I)+N233* (RPHI 


AE(1,2)=C21+C22+PN4R 
AE(1,3)=0.0D0 


I)**2) 
(I)**2) 


C41=N3+N113*(V(I) **2)+2. ODO*N133*V (I) *RPHI (I) +N223* (R(I)**2) 
C42=2 .0D0*N233*R(I)*RPHI(I) 


AE(1,4)=C41+C42 


D11=Y1+3.0D0*Y111*(V(I) **2)+¥122* (R(I) **2)+2.0D0*Y112*V(I)*R(I) 
D12=2.0D0*Y113*V(1I) *RPHI(I)+Y133* (RPHI (I) **2) 


AE(2,1)=D11+D12+PY4V 


D21=2.0D0*Y122*V(I)*R(1I)+(Y2-M*UND) +3 .0D0*Y222* (R(I)**2) 
D22=2.0D0*Y223*R (I) *RPHI(1)+Y233* (RPHI (I) **2)+¥112*(V(I)**2) 


AE (2,2) =D21+D22+PY4R 
AE(2,3)=0.0D0 
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D41=Y3+Y113*(V(I) **2)+2.0D0*Y133*V(I) *RPHI(1)+¥223* (R(1I) **2) 
D42=2.0D0*Y233*R(1I) *RPHI (I) 


AE (2,4) =D41+D42 


E11=K1+3.0D0*K111*(V(I) **2)+K122* (R(I)**2)+2.0D0*K112*V(I)*R(I) 
E12=2.0D0*K113*V(I) *RPHI(1I)+K133* (RPHI(I)**2) 


AE(3,1)=E11+E12+PK4V 


E21=2.0D0*K122*V(1)*R(1)+(K2+M*ZG*UND) +3 .0D0*K222* (R(I) **2) 
E22=2.,.0D0*K223*R(1I) *RPHI (1) +K233* (RPHI (I) **2) +K112* (V(1) **2) 


AE(3,2)=E21+E22+PK4R 
AE(3,3)=0.0D0 


DGZ1=GM*DCOS(RPHI(1I)) 
DGZ2=BM* (DTAN (RPHI (I) 


1) )+0.5D0*BM* (DTAN(RPHI 
PB Teg ( 
DGZ=DGZ1+DGZ2 
) 
2 
* 


I) ))**2*DCOS(RPHI ( 
((1.0D0/DCOS(RPHI(I)))* 


*2) *DSIN(RPHI ( 


eS 


( )) 

I ))) 

DGZ=GM*DCOS(RPHI(I 

ene yi2) 72) 
I 


) 
)+2.0D0*K133*V (I) *RPHI (1) +K223* (R(I)**2) 
E42=2.0D0*K233*R(I) *RP 


HI(1) -WG*DGZ 


AE (3,4) =E41+E42 
AE(4,1)=0.0D0 
AE(4,2)=0.0D0 
AE(4,3)=1.0D0 
AE(4,4)=0.0D0 


This 1s a subroutine that utilizes the EISPACK LIBRARY [Ref. 20] 
Purpose: Calls the recommended sequence of subroutines from the 


Eigensystem subroutine package (EISPACK) to find the 
eigenvalues ou nue real generalized eigenproblem: 


= (LAMBDA) *BX 
Arguments: 

NM -Row dimension of the two-dimensional array parameters as 
declared in the calling program Dimension statement 
(Input) 

N -The order of the matrices A and B. (Input) 

A -Contains a real general matrix (Input) 

B -Contains a real general matrix (Input) 


MATZ -Path indicator. If= 0, eigenvalues are found. (Input) 

ALFR -Real parts of the numerators of the eigenvalues. (Output) 

ALFI -Imaginary parts of the numerators of the eigenvalues. 
(Output) 

BETA -Contains the denominators of the eigenvalues, which are 
given by the ratios: (ALFR+I*ALFI)/BETA. (Out put ) 

Z -For eigenvector usage if MATZ not equal to suc (Output) 

IERR -Error completion code. (Output) 


OAleRGG( 4. 4 ,AE,BE,ALFR,ALFI,BETA,0,Z, IERR) 
DO 600 J=1, 
IF (BETA(J). “NE.0.0D0) THEN 
EVALR(J)=ALFR(J) /BETA(J) 
EVALI (J) =ALFI (J) /BETA(J) 
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ELSE 
EVALR (J) =DMACH (2) 
EVALI (J) =DMACH (2) 
ENDIF 
600 CONTINUE 
WRITE (25,25) EVALR(1) ,EVALI(1) ,EVALR(2) , EVALI(2) 
WRITE (30,25) EVALR(3) ,EVALI(3) , EVALR(4) , EVALI(4) 


* Calculate the damping ratio and the natural frequency 
* 


Z1=-EVALR(1) /(DSQRT(EVALI (1) **2+EVALR (1) **2) ) 
Z2=-EVALR(2) / (DSQRT(EVALI (2) **2+EVALR (2) **2) ) 
Z3=-EVALR(3) / (DSQRT(EVALI (3) **2+EVALR (3) **2) ) 
Z4=-EVALR(4) / (DSQRT(EVALI (4) **2+EVALR(4) **2) ) 
WN1=DSQRT(EVALI (1) **2+EVALR (1) **2) 
WN2=DSQRT (EVALI (2) **2+EVALR (2) **2) 
WN3=DSQRT (EVALI (3) **2+EVALR (3) **2) 
WN4=DSQRT (EVALI (4) **2+EVALR (4) **2) 
WRITE(35,25)Z1,Z2,Z3,24 
WRITE (36, 25)WN1,WN2,WN3,WN4 
25 FORMAT(4E15.5) 
500 CONTINUE 
RETURN 
END 


SUBROUTINE EIGTH2(VDEL,V,R) 


INCLUDE ’HYDCOMN. for’ 
INTEGER I,J,IERR 
REAL*8 AE(2,2) ,BE(2,2) ,BETA(2) , VDEL(300) ,V(300) ,R(300) 
REAL*8 EVALR(2) ,EVALI(2) ,EVAL(2) ,ALFR(2) ,ALFI(2) ,Z2(2,2) 
EXTERNAL DMACH,RGG 
DO 500 I=1,imax 

* SET UP THE B MATRIX AND RESET EACH I 
UP=UND* (WP+TAU* ((V(I)+XP*R(I))**2+CPV*V(I)+CPR*R(I) )) 
V9=V(1I)*0.99D0 
R9=R(I)*0.99D0 
UPV9=UND* (WP+TAU* ( (V9+XP*R(I) ) **2+CPV*V9+CPR*R(1I) ) ) 
UPRQ=UND* (WP+TAU* ( (V(1)+XP*RQ) **2+CPV*V(I)+CPR*RQ) ) 
JP=(U*UP) / (NP*D) 
JPV9=(U*UPV9) / (NP*D) 
JPR9=(U*UPRQ) / (NP*D) 
KT=KTA+ (KTB*JP) 
KTV9=KTA+(KTB*JPV9 ) 
KTRO=KTA+(KTB*JPRQ) 
UR=UP*EPSILON*((1.0D0+8.ODO*K*KT) /(PI*(JP**2) ))**0.5 


URV9=UPV9*EPSILON* ( (1.0D0+8.0D0*K*KTVQ) / (PI* (UPV9**2) ))**0.5 
URR9=UPR9*EPSILON* ((1.0D0+8.O0D0*K*KTRQ) /(PI* (UPR9**2) ))**0.5 
VR=GAMMA*V (I) +CDR*R(1I)+CDRRR*(R(1) **3)+CDRRV*V(I) *(R(1I)**2) 
VRV9=GAMMA*V9+4CDR*R (1) +CDRRR* (R(1) **3)+CDRRV*V9* (R(1) **2) 


VRR9=GAMMA*V (1) +CDR*R9+CDRRR* (R9**3) +CDRRV*V (I) * (RO**2) 
VD=(UR**2+VR**2) **0.5 
VDV9=(URV9**2+VRV9**2) **0 5 
VDR9=(URRO**2+VRR9**2) **0.5 
ALPHAR=VDEL (1) +DATAN(VR/UR) 
ALPHARV=VDEL (1) +DATAN (VRV9/URV9 ) 
ALPHARR=VDEL (1) +DATAN(VRR9/URRQ) 
FLAM=(6.13D0*LAMBDA) / (LAMBDA+2.25D0) 

FNN=FLAM* (AR/ (L**2) )*( (VD) **2) *DSIN(ALPHAR) 
FNNV9=FLAM* (AR/ (L**2) ) *((VDV9) **2) *DSIN(ALPHARV) 
FNNRO=FLAM* (AR/(L**2) ) *((VDRQ) **2) *DSIN(ALPHARR ) 
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% 


Y4=-(1.0D0+AH) *FNN*DCOS(VDEL (I) ) 
Y4V9=- (1.0D0+AH) *FNNV9*DCOS (VDEL 
Y4R9=- (1.0D0+AH) *FNNR9*DCOS(VDEL 
PY4V=(Y4-Y4V9) /(V(I) -V9) 
PY4R=(Y4-Y4R9)/(R(1) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS(VDEL (1) ) 
N4V9=- (XR+AH*XH) *FNNV9*DCOS(VDEL (1) ) 
N4R9=- (XR+AH*XH) *FNNR9*DCOS(VDEL (1) ) 
PN4V=(N4-N4V9) /(V(I) -V9) 
PN4R=(N4-N4R9)/(R(1)-R9) 


ny, gm, 


WN= ( (WG*GM) / (0.000021D0) )**0.5 
Equation 4.4b 


BE(1,1)=0.00035245D0 
BE(1,2)=0.000875D0 
BE(2,1)=0.014969D0 
BE (2,2) =0.00035245D0 


Equation 4.4c-g 
C11=N14+3*N111* (V(1) **2)+N122* (R( 1) **2)+2.0D0*N112*V(1)*R(T) 
AE(1,1)=C11+PN4V 


C21=2 .0D0*N122*V(I) *R(1)+(N2-M*XG*UND) +3 .0D0*N222*(R(I)**2) 
C22=N112*(V(I) **2)+PN4R 


AE(1,2)=C21+C22 
D11=Y1+3.0D0*Y111*(V(I)**2)+Y¥122*(R(1)**2)+2.0D0*Y112*V (1) *R(1) 
AE (2,1)=D11+PY4V 


D21=2.0D0*Y122*V(I)*R(1)+(Y2-M*UND) +3 .0D0*Y222* (R(I) **2) 
D22=Y112* (V(I)**2)+PY4R 


AE (2,2) =D21+D22 


CALL RGG(2,2,AE,BE,ALFR,ALFI,BETA,0,Z, IERR) 
DO 600 J=1,2 
IF (BETA(J) .NE.0.0D0O) THEN 
EVALR(J)=ALFR(J) /BETA(J) 
EVALI (J) =ALFI(J) /BETA(J) 
ELSE 
EVALR (J) =DMACH (2) 
EVAL I (J) =DMACH (2) 
ENDIF 
600 CONTINUE 
WRITE (26,25) EVALR(1) ,EVALI(1) ,EVALR(2) , EVALI(2) 
25 FORMAT(4E15.5) 
500 CONTINUE 
RETURN 
END 


SUBROUTINE EIGTH3(VDEL,RPHI) 


INCLUDE ’HYDCOMN. for’ 
INTEGER I,J,IERR 
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* 


ga 0 


REAL*8 AE(2,2) ,BE(2,2) ,BETA(2) , VDEL(300) , RPHI(300) 
REAL*8 P(300) ,EVALR(2) ,EVALI(2) ,EVAL(2) ,ALFR(2) ,ALFI(2) ,Z(2,2) 


EXTERNAL DMACH,RGG 
DO 500 I=1,imax 
P(I)=0.0D0 


SET UP THE B MATRIX AND RESET EACH I 


600 


25 
500 


WN=( (WG*GM) / (0 .000021D0) )**0.5 
IF(FN.LE.O.1) THEN 

KAPPA=0.1D0 

ELSE IF (EN. GE.0.2) THEN 

KAPPA=0 .2D0 

ELSE 

KAPPA=FN 

ENDIF 
KPHIDOT=-0.000021D0*KAPPA*WN 


Equation 4.7b 


BE(1,1)=0.000021D0 
BE(1,2)=-KPHIDOT 
BE(2,1)=0.0d0 
BE(2,2)=1.0D0 


Equation 4.7c-e 


AE(1,1)=0.0d0 
DGZ1=GM*DCOS (RPHI (I) )+0.5D0*BM* (DTAN(RPHI 
DGZ2=BM* (DTAN(RPHI (I) )*((1.0D0/DCOS(RPHI ( 
DGZ=DGZ1+DGZ2 
DGZ=GM*DCOS (RPHI(1) ) 

AE(1,2)=K3-WG*DGZ 

AE(2,1)=1.0D0 

AE(2,2)=0.0D0 

CALL RGG(2, 2,AE,BE,ALFR,ALFI,BETA,0,Z, JERR) 

DO 600 J=1,2 
IF (BETA(J) .NE.0.ODO) THEN 

EVALR(J)=ALFR(J) /BETA(J) 
EVALI(J)=ALFI (J) /BETA(J) 
ELSE 
EVALR (J) =DMACH (2) 
EVALI (J) =DMACH (2) 
ENDIF 

CONTINUE 

WRITE (27,25) EVALR(1) , EVALI(1) ,EVALR(2) , EVALI(2) 
Z1=-EVALR(1)/(DSQRT(EVALI(1) **2+EVALR(1) **2) ) 
Z2=-EVALR (2) / (DSQRT(EVALI (2) **2+EVALR (2) **2) ) 
WN1=DSQRT(EVALI (1) **2+EVALR(1) **2) 
WN2=DSQRT (EVALI (2) **2+EVALR(2) **2) 

WRITE (37, 25)Z1,Z2,WN1,WN2 

FORMAT (4E15.5) 

CONTINUE 

RETURN 

END 


(I) )) 
T)))* 
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APPENDIX B 
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Fn = 0.3m and a = 0.6 
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B.2: Effect of GM on Coupled Model Roll Eigenvalues; Imaginary Part 
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B.3: Effect of GM on Coupled Model Upper Steering Eigenvalues 
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GM = 0.3m and a@a = 0.6 
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B.4: Effect of GM on Coupled Model Lower Steering Eigenvalues 
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GM = 0.3m and a = 0.6 
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B.5: Effect of Fn on Decoupled Model Damping Ratio 
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B.6: Effect of Fn on Coupled Model Roll Eigenvalues 
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B.7: Effect of Fn on Coupled Model Upper Steering Eigenvalues 
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B.8: Effect of Fn on Coupled Model Lower Steering Eigenvalues 
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B.9: Effect of Fn on Coupled Model Roll Eigenvalues 
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B.10: Effect of Fn on Coupled Model Root Locus 
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GM = 0.3 m and Fn = 0.3 
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B.11: Effect of Speed Loss on Coupled Model Roll Eigenvalues 
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B.12: Effect of Speed Loss on Coupled Model Root Locus 
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B.13: Effect of Speed Loss on Coupled Model 
Eigenvalues 
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B.14: Effect of Speed Loss on Coupled Model Lower Steering 
Eigenvalues 
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B.15: Effect of Speed Loss on Coupled Model Sway Velocity with GM = 1.0m 
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B.16: Effect of Speed Loss on Coupled Model Sway Velocity with GM = 3.0 m 
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B.17: Effect of Speed Loss on Coupled Model Roll Angle with GM = 1.0 m 
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B.18: Effect of Speed Loss on Coupled Model Roll Angle with GM = 3.0 m 
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B.19: Effect of Speed Loss on Coupled Model Root Locus with GM = 1.0m 
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B.20: Effect of Speed Loss on Coupled Model Root Locus with GM = 3.0 m 
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